/* José Filipe M Santos - 7277502 */

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <gmp.h>
#include <math.h>
#include<string.h>
#include<pthread.h>
#include<unistd.h>
#include<time.h>

#define DIGITOS 10000000

/* algoritmo de borwein com convergencia quadratica */

int i = 0;										/* numero de iteracoes 						*/
unsigned long int um = 1, dois = 2;							/* constantes utilizadas 					*/
pthread_t tid[5];									/* ids para as threads 						*/
mpf_t a0, b0, p0, sqrta0, intermedio1, intermedio2, pi, diff, zero;			/* a0, b0, p0: variaveis da formula 				*/
											/* srtqa0: raiz de a0     					*/
											/* intermedio i: valores intermediarios utilizados nos calculos	*/
											/* pi, diff, zero: valores utilizados para verificar cvergencia	*/

pthread_mutex_t lockfile;								/* mutex para sincronizar impressao				*/

FILE *ifp;										/* arquivo para impressao dos dados				*/

void* tb00(void *arg){									/* threads para os calculos relativos a b0 */

	mpf_add_ui(intermedio1, b0, um);	
	mpf_mul(intermedio1, intermedio1, sqrta0);
	return NULL;

}

void* tb01(void *arg){

	mpf_add(intermedio2, a0, b0);
	return NULL;

}

void* ta00(void *arg){									/* threads para os calculos relativos a a0 */

	mpf_add_ui(intermedio1, a0, um);
	return NULL;

}

void* ta01(void *arg){

	mpf_mul_ui(intermedio2, sqrta0, dois);
	return NULL;

}

void* tp00(void *arg){									/* threads para os calculos relativos a p0 */

	mpf_add_ui(intermedio1, a0, 1);
	return NULL;

}

void* tp01(void *arg){

	mpf_add_ui(intermedio2, b0, 1);
	return NULL;

}

void* tprint(void *arg){

	pthread_mutex_lock(&lockfile);

	mpf_out_str (ifp, 10, DIGITOS, p0);					/* inserir pi a cada iteracao no arquivo */
	fprintf(ifp, " \n");

	pthread_mutex_unlock(&lockfile);

	return NULL;

}

void chama_tb(void){									/* função que paraleliza os calculos de b0 */


	pthread_create(&(tid[0]), NULL, &tb00, NULL);
	pthread_create(&(tid[1]), NULL, &tb01, NULL);

	pthread_join(tid[0], NULL);
	pthread_join(tid[1], NULL);
}

void chama_ta(void){									/* função que paraleliza os calculos de a0 */


	pthread_create(&(tid[0]), NULL, &ta00, NULL);
	pthread_create(&(tid[1]), NULL, &ta01, NULL);

	pthread_join(tid[0], NULL);
	pthread_join(tid[1], NULL);
}

void chama_tp(void){									/* função que paraleliza os calculos de p0 */


	pthread_create(&(tid[0]), NULL, &tp00, NULL);
	pthread_create(&(tid[1]), NULL, &tp01, NULL);

	pthread_join(tid[0], NULL);
	pthread_join(tid[1], NULL);
}

void chama_print(void){

	pthread_create(&(tid[0]), NULL, &tprint, NULL);

	pthread_join(tid[0], NULL);

}

void borwein(){										/* algoritmo de borwein */

	mpf_set_default_prec(DIGITOS * sizeof(int));					/* estabeler precisao padrao */
	mpf_init(a0);									/* inicializar variaveis */
	mpf_init(sqrta0);
	mpf_init(b0);
	mpf_init(p0);
	mpf_init(pi);
	mpf_init(diff);
	mpf_init(intermedio1);
	mpf_init(intermedio2);
	mpf_init(zero);

	mpf_sqrt_ui(a0, dois);								/* fixando valores iniciais (a0 = sqrt(2) e p0 = 2 + sqrt(2))*/
	mpf_sqrt_ui(p0, dois);
	mpf_add_ui(p0, p0, dois);

	ifp = fopen("pi.txt", "w+");						/* arquivo para escrever os dados */

	do{

		mpf_set(pi, p0);							/* pi = p0 atual */
		mpf_sqrt(sqrta0, a0);

		chama_tb();

		mpf_div(b0, intermedio1, intermedio2); 					/* b[n+1] = (sqrt(a[n]) * (1 + b[n])) / (a[n] + bn[n]) */

		chama_ta();

		mpf_div(a0, intermedio1, intermedio2);					/* a[n+1] = a[n] + 1 / 2 x sqrt(a[n]) */

		chama_tp();

		mpf_div(intermedio2, intermedio1, intermedio2);
		mpf_mul(intermedio1, b0, p0);
		mpf_mul(p0, intermedio1, intermedio2);					/* p[n+1] = (p[n] * b[n+1])*((1 + a[n+1])/(1 + b[n+1])) */

		mpf_sub(diff, p0, pi);							/* tirar a diferenca, se for 0, está pronto */
		chama_print();

		++i;

		


	}while(mpf_cmp(zero, diff) != 0);						/* iterar ate diff == 0 */

	fprintf(ifp, " iteracoes: %d\n", i);						/* inserir numero total de iteracoes */
	fclose(ifp);	
	return;

}

int main(){

	time_t inicio, fim;

	time(&inicio);
	pthread_mutex_init(&lockfile, NULL);
	borwein();
	time(&fim);
	printf("RELATORIO DE EXECUCAO:");
	printf("%d ITERACOES EM %lfs\n", i, difftime(fim, inicio));

	return 0;
}
